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A propagation method for timc-dcpendcnt Schrodingcr equations with an exphcitly time-dependent Hamil- 
tonian is developed where time ordering is achieved iteratively. The exphcit time-dependence of the time- 
dependent Schrodinger equation is rewritten as an inhomogeneous term. At each step of the iteration, the 
resulting inhomogeneous Schrodinger equation is solved with the Chebychev propagation scheme presented 
in J. Chem. Phys. 130, 124108 (2009). The iteratively time-ordering Chebychev propagator is shown to be 
robust, efficient and accurate and compares very favorably to all other available propagation schemes. 



I. INTRODUCTION 

The dynamics of the interaction of matter with a 
strong radiation field is described by time-dependent 
Schrodinger equations (TDSEs) where the Hamiltonian 
is explicitly time-dependent. This description is at the 
core of the theory of harmonic generation^ pump- 
probe spectroscopy,— and coherent controL^^ Typically, 
an atom or molecule couples to a laser pulse via a dipole 
transition, 



The difficulty of simulating explicitly time-dependent 
Hamiltonians, emerges from the fact that the commuta- 
tor of the Hamiltonian with itself at different times does 
not vanish 1^ 



[H(ii),H(t2)]- ^0. 



(2) 



Formally, this effect is taken into account by time order- 
ing such that the time evolution is given by 



U(r,0) ^Te^^^o^H 



(t) dt 



(3) 



H(t) = Ho + ^(i)A, 



(1) 



with E{t) the time-dependent electromagnetic field, caus- 
ing the explicit time-dependence of the Hamiltonian. 
Simulating these light-matter processes from first princi- 
ples imposes a numerical challenge. Realistic simulations 
require efficient procedures with very high accuracy. 

For example, in coherent control processes, interaction 
of quantum matter with laser light leads to constructive 
interference in some desired channel and destructive in- 
terference in all other channels. In time-domain coherent 
control such as pump-probe spectroscopy, wave packets 
created by radiation at an early time interfere with wave 
packets generated at a later time. This means that the 
relative phase between different partial wave packets has 
to be maintained for long time with high accuracy. As 
a result, numerical methods designed to simulate such 
phenomena have to be highly accurate, minimizing the 
errors in both amplitude and phase. 
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The effect of time ordering is to incorporate higher order 
commutators into the propagator U(T, 0). For strong 
fields E{t) and fast time-dependences the convergence 
with respect to ordering is slow. Methods to incorporate 
the second order Magnus termi have been developed ei- 
ther in a low order polynomial expansion^iS or as a split 
exponential 

A quantum dynamical propagator that fully accounts 
for time ordering is given by the {t,t') method-ii It 
is based on rewriting the Hamiltonian in an extended 
Hilbert space where an auxiliary coordinate, t' , is added 
and terms such as E{t')fi are treated as a potential in this 
degree of freedom. The Hamiltonian thus looses its ex- 
plicit dependence on time t, and can be propagated with 
one of the available highly accurate methods for solving 
the TDSE with time-independent Hamiltonian^^ 

Most of the vast literature on the interaction of mat- 
ter with time-dependent fields in genera&i^rJ^ and on 
coherent control in particulari^ii^"— ignores the effect of 
time ordering. Popular approaches include Runge-Kutta 
schemes,— 1^°'^^ the standard Chebychev propagator with 
very small time step^^ and the split propagator j ^"^'^^:^^ 
Naively it is assumed that if the time step is small enough 
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the calculation with an explicit time-dependent Hamilto- 
nian can be made to converge. The diiEculty is that this 
convergence is very slow - second order in the time step 
if the Hamiltonian is stationary in the time interval and 
third order if the second order Magnus approximation is 
usedi^i^ Additionally in many cases the error accumu- 
lates in phase^i^ so that common indicators of error such 
as deviation from unitarity are misleading. 

In order to obtain high quality simulations of explic- 
itly time dependent problems a new approach has to be 
developed. The ultimate (i, t') method cannot be used in 
practice since it becomes prohibitively expensive in real- 
istic simulations. On the other hand we want to maintain 
the exponential convergence property of spectral decom- 
position such as the Chebychev propagator. The solution 
is an iterative implementation of the Chebychev propa- 
gator for inhomogcneous equations such that it can over- 
come the time ordering issue. 

The paper is organized as follows. The formal solution 
to the problem is introduced in Section [TTl The TDSE 
for an explicitly time-dependent Hamiltonian is rewrit- 
ten as an inhomogeneous TDSE. The inhomogeneity is 
calculated itcratively and converges in the limit of many 
iterations. At each step of the iteration, an inhomogene- 
nous TDSE is solved by a Chebychev propagator which is 
based on a polynomial expansion of the inhomogeneous 
term^S. The resulting algorithm is outlined explicitly in 
Section IIIII and applied to three different examples in 
Section HVl Its high accuracy is demonstrated and its ef- 
ficiency is discussed in comparison to other approaches. 
Section |V] concludes. 



II. FORMAL SOLUTION 

The Hamiltonian, H, describing the interaction of a 
quantum system with a time-dependent external field 
topically consists of a field-free, time-independent part. 
Ho, and an interaction term, W(t) = piE{t). The TDSE 
for such a Hamiltonian (setting /i = 1), 

^^^\m) = {Ho + \N{t))m))^ (4) 

is solved numerically by dividing the overall propagation 
time [0,r] into short time intervals [i„,t„^_i], each of 
length At. A two-stage approach is employed. First, 
the formal solution of the TDSE is considered. The term 
arising from the explicit time-dependence of the Hamilto- 
nian is approximated iteratively. The iterative loop thus 
takes care of the time ordering. Second, at each step 
of the iteration, an inhomogeneous Schrodinger equation 
is obtained. It is solved with the recently introduced 
Chebychev propagator for inhomogeneous Schrodinger 
equations 1^ 



A. Iterative time ordering 

The TDSE, Eq. (|4]), is rewritten to capture the time- 
dcpendcncc within the interval [tn,tn+i\, 

= (Ho + W„)lV^W> + (W(t)-W„)|^(i)> . (5) 

Here, W„ is the value of W(t) at the midpoint of the 
propagation interval, W„ = \N j , xhe formal 

solution of Eq. ^ is given by 

|^(i))=e-^""(*-*")K.(f„)>- 

ij^ e-*""(*--)V„(T)|V^(r))dr, (6) 

where H„ = Hq + W„ denotes the part that is indepen- 
dent of time in and V„(i) = W(i) — Wn the 
time-dependent part. Eq. ([6]) is subjected to an iterative 
loop, 

*^*e-^""^*-^)V„(r)|^fc_i(r))dr, (7) 

The solution at the fcth step of the iteration, \ipk), is 
calculated from the formal solution, Eq. ©, by replacing 
IV'fc) in the second term on the right-hand side of Eq. ([6]) 
by \ipk_i) which is known from the previous step. 

In this approach, time ordering is achieved by converg- 
ing \ipk-i) to \ipk) as the iterative scheme proceeds. This 
is equivalent to the derivation of the Dyson series. Start- 
ing from the equation of motion for the time evolution 
operator, 

j^O(t,0) = H(t)U(t,0), 
the formal solution for the time evolution operator, 

U(t,0) ^-2 f H(ti)U(ti,0)dii, (8) 
Jo 

is iteratively inserted in the right-hand side, i.e. 

U(t,0)--z / / ' H(ii)H(i2)U(t2,0)dt2dii, 

Jo Jo 

U(t,0) = -z / / ... / 

Jo Jo Jo 

H{ti)f\{h) . . . H(t„)U(t„, 0)dtn . . . dt2dti , 

where U(t„,0) goes to 1 as t„ becomes smaller and 
smaller. Our formal solution, Eq. ([6]) is equivalent to 
Eq. ([S]). An alternative approach to time ordering is 
given by the Magnus expansion which is based on the 
group properties of unitary time evolutiouii In the limit 
of convergence, the Magnus and the Dyson series are 
completely equivalent, but low-order approximations of 
the two differ Our iterative scheme corresponds to the 
limit of convergence (with respect to machine precision) . 
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B. Equivalence to an inhomogeneous TDSE 

Differentiating Eq. ([7]) with respect to time, an inlio- 
mogcneous Sclirodinger equation at each step k of the 
iteration is obtained, 



dt 



The inhomogencity is given by 



(9) 



(10) 



Eq. ([9]) can be solved by approximating the inhomoge- 
neous term globally within [tn,tn+i], i-C. by expanding 
it into Chebychev polynomials, 

m— 1 

^P,©|$fc_l,,). (11) 

Pfc-i.j denotes the Chebychev polynomial of order j with 
expansion coefficient \^k-i,j), and t = 2{t — tn)/^t — 1 
with t G [<„,i„+i] is a rescaled tme.— 

The expansion coefficients, in Eq- (fTT|) are 

given by 



1 VT^ 



(12) 



Since |$fe-i(i)) is known at each point in the interval 
and in particular at the zeros, ti, of the mth Chebychev 
polynomial, the integral in Eq. (|12p can be rewritten by 
applying a Gaussian quadrature;^ yielding 



2 - Sjo 



m — 1 



(13) 



i=0 



Due to the fact that the Chebychev polynomials can be 
expressed in terms of cosines, Eq. (jl3p is equivalent to a 
cosine transformation. Thus the expansion coefficients, 
|$fc_ij), can easily be obtained numerically by fast co- 
sine transformation. 

The expansion into Chebychev polynomials, if con- 
verged, is equivalent to the following alternative expan- 
sion. 



m — 1 



j'=0 



(14) 



Once the coefficients of the Chebychev expansion, 
|^fe-i.j)j are known, the transformation described in Ap- 

pendix 1^ is used to generate the coefficients |$)^:'_]^) in 
Eq. HID. 

Approximating the inhomogeneous term by the right- 
hand side of Eq. , the formal solution of Eq. (O can 
be writtenS^ 

\Mt)) = E ^^-r^l^^i) + - (15) 

j=o 



where the lA^"*^!) are obtained recursively, 

\\tll)=mn)). (16) 

1 < j < m . 
Fm is a function of H„ and is given by 
F,„= (17) 

{-lHn{t-tn)y 




Taking the derivative of Eq. (jTS]) with respect to time, the 
inhomogeneous Schrodinger equation is recovered after 
some algebra.— 

Alternatively, Eq. pn]) can be inserted into Eq. ([7]), 
replacing by its polynomial approximation, 

Eq. ini), 

\m) = er'^-'\m) + 

rn-l „t j 



3=0 



J! 



(without any loss of generality, t„ has been set to zero). 
Defining 



r 



(19) 



and integrating Eq. p9|) by parts, one obtains 



= (-^H„)-^ ( - |a ) , (20) 



I < j < m — 1 , 
By induction, it follows that 



(21) 



«, = (-^H„)-(^■+^' ( e-^^"* - E ) • (22) 



a=0 



Defining 



(23) 



a=0 



Eq. pS]) becomes 



|V>(t)) = e-""*|^(0)) + E F.+il*^''^> , (24) 

3=0 



which was shown to be equivalent to Eq. p^ .— 

The algorithm for solving the TDSE with explicitly 
time-dependent Hamiltonian is thus based on evaluating 
the integral of the formal solution, Eq. ([S]), in an iterative 
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fashion. At each step k of the iteration, the inhomoge- 
neous TDSE, Eg. is solved by applying the propaga- 
tor of Ref. lla within each short time interval [in, ^n+i]- 

Once convergence with respect to the iteration k is 
reached, the inhomogeneous term becomes constant with 
respect to k. 



III. OUTLINE OF THE ALGORITHM 

We assume that the action of the Hamiltonian on a 
wavefunction can be efficiently computedJ^ Then the 
complete propagation time interval [0, T] is split into 
small time intervals, [t„,tri+i]- For each time step 
[tmtn^i], the implementation of the Chcbychcv propa- 
gator with iterative time ordering involves an outer loop 
over the iterative steps k for time ordering and an in- 
ner loop over j for the solution of the (inhomogeneous) 
Schrodinger equation for each k. 

1. Preparation: Set a local time grid {r;} for each 
short-time interval [tmtn+i]. In order to calcu- 
late the expansion coefficients of the inhomoge- 
neous term by cosine transformation, the Nt sam- 
pling points {ti} are chosen to be the roots of the 
Chebychcv polynomial P/v^ of order N^. The num- 
ber of sampling points, Nt, is not known in advance. 
One thus has to provide an initial guess and check 
below, in step 3.i, that it is equal to or larger than 
the number of Chebychev polynomials required to 
expand the inhomogeneous term. 



Nt> m. 



(25) 



If Nt is much larger than m, it is worth to decrease 
it (subject to the bound of Eq. and to recal- 

culate the {ti}. The number of propagation steps 
within [t„,t„+i] is then reduced to its mininum. 

2. The propagation for fc = solves the Schrodinger 
equation for the time-independent Hamiltonian 
H„ = Ho+W„, 

'^\Mr)) = (Ho + W„)IV^(t)), 

with initial condition \4'o{t = t„)) ~ |i/'(t„)). A 
standard Chebychev propagator is employed to this 
end. Note that for k = 0, the same time grid {t;} 
needs to be used as for fc > because the inho- 
mogeneous term for fc = 1 is calculated from the 
zeroth order solution, \'ipo{t)). Since the {r;} are 
not equidistant, the Chebychev expansion coeffi- 
cients of the standard propagator, g"*'^"'^'^', need 
to be calculated for each time step within [t„, 
where At; = t;+i - t;, / = 1, Nt - 1. 

3. The fc > propagation solves an inhomogeneous 
Schrodinger equation, cf. Eq. ([5]), with the ini- 
tial condition \iph{t = t„)) = \i^o{t = t„)) = 



\ip{tn))- This is achieved by the Chebychev propa- 
gator for inhomogeneous Schrodinger equations 
i.e. Eq. (fT5|l . and involves the following steps: 

(i) Evaluate the inhomogeneous term, 

|$fc-l(T)) - -z(W(t) - W„)|Vfc-l(T)). 

(ii) Calculate the expansion coefficients of the 
inhomogeneous term, cf. Eqs. (jll[) and 

The Chebychev expansion coefficients 
|4>fc_ij) are obtained by cosine transforma- 
tion of |$fc_i(r))^ The coefficients |$^^2i) 
are evaluated from the Chebychev expansion 
coefficients \^k~i,j) using the recursive rela- 
tion given in Eqs. (jA14[) and (jAlSp . The order 
m of the expansion is chosen such that ratio 
of the smallest to the largest Chebychev coeffi- 
cient becomes smaller than the specified error 



1$ 



fc-l,?ri+l| 



< e. 



fe-l,0| 



(26) 



To obtain high accuracy, e may correspond to 
the machine precision.— 

(iii) Calculate the Chebychev expansion coeffi- 
cients of F„i, cf. Eq. ([TT]) . also by cosine 
transformation. The number of terms in 
this Chebychev expansion is also determined 
by the relative magnitude of the coefficients, 
analogously to Eq. (|26|) . 

(iv) Determine all required in Eq. (fTSj) by 
evaluating Eq. (fTHll . 

(v) Construct the solution |V'fc(t = tn+i)) accord- 
ing to Eq. m 



4. Convergence is reached when \ipk-i(tn+i)) and 
I'ijjkitn+i)) become indistinguishable, 

||^.'=-i(i„+i)-|V''=(<„+i)|| <e, 

and the desired solution of the Schrodinger equa- 
tion with explicitly time-dependent Hamiltonian is 
obtained, |V'(t„+i)) = \ipkitn+i)) ■ 

The only parameter of the algorithm is the pre- 
specified error e. It determines the number of iterative 
terms fc and the order of the inhomogeneous propagator 
m. Furthermore, to execute the algorithm, the user has 
to provide, besides e, an initial guess for the number of 
sampling points of the local time grid, Nt. 



IV. EXAMPLES 

We test the accuracy and efficiency of the algorithm for 
three examples of increasing complexity. The first two 
examples, a driven two-level atom and a linearly driven 
harmonic oscillator, are analytically solvable. We can 
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therefore compare the numerical to the analytical solu- 
tion and establish the accuracy of the Chebychev propa- 
gator with iterative time ordering. For the third example, 
wave packet interferometry in two oscillators coupled by 
a field, no analytical solution is known. The Chebychev 
propagator with iterative time ordering thus serves as a 
reference solution to which less accurate methods can be 
compared. 



A. Driven two-level atom 

The Hamiltonian for a two-level atom driven reso- 
nantly by a laser field in the rotating- wave approximation 
reads^^ 



"-Ui^W j 



where the field is of the form 



E{t) = -EoS{t) , 



(27) 



(28) 



and S{t) denotes the envelope of the field. We take the 
strength of the transition dipole to be /i = f a.u., the final 
propagation time T = 9000 a.u., and the shape function 



S{t) = sin^ 



Analytically, the time evolution of the amplitudes. 



is obtained as 



c^"'^(<) = cos 



'it) 



1 { T /2TTt 

I ^ f T . /2TTt 

— uEn t sm 

4^ " V 27r \ T 



(29) 

(30) 

(31) 
(32) 



Initially the two-level system is assumed to be in the 
ground state, Cg(< = 0) = 1, Ce(t = 0) = 0. The pulse 
amplitude is chosen to yield a 7r-pulse, such that c^^{t = 
T) = 0, c|"^(t = T) = 1. 

Defining at each time step the errors. 



and 



eUt)^\\cr{i)?-\cM^\ : 



,{t)^\i-m)m))\ 



(33) 



(34) 



we measure the deviation of the numerical from the an- 
alytical solution and the deviation of the norm of \i^{t)) 
from unity. The time evolution of esoi(i) and enorm(^) 
is shown in Fig. [1] for the Chebychev propagator with 
iterative time ordering. The maximum errors occuring 
during the propagation, e 



max 
norm 



and 



■sol 



marized in Table[Tl For time steps, At = f 



are also sum- 
„+i -t„, up to 



1x10"° 
IxlO"' 
ixio'" 



1x10 
1x10"' 
1x10"' 
1x10' 



1x10 



1x10 
0? 1x10 
1x10 
1x10 



1 1 ' 1 r 


1 1 1 1 1 1 1 1 1 1 r 


1 ' 


r *Sc « 


^_ • * — - 






><-«At= 40 a.u. \ 






At = 300 a.u. \ 




? / . , 1 


At = 900 a.u. V. 
, . , 1 , . , 1 , . , 






4000 6000 
time t [a.u.] 



FIG. 1. (color online) Strongly driven two-level atom prop- 
agated with iteratively time ordering Chebychev propagator: 
error of the solution, eaoi(i) (a) and deviation of the norm of 
\ii{t)) from unity, enorm(t) (b). 



about T/lOO, the maximum errors occuring during the 
propagation, ef^^, are of the order of 10~^^. If the time 
step is further increased to about T/10, the maximum 

is 



sol 



errors are of the order of 10"'^. The increase in e; 
accompanied by an increase in ej^jj^^ as the time steps 
become larger, cf. Table HI The deviation from unitarity 
indicates that the error is due to the Chebychev expan- 
sion of the time evolution which becomes unitary only 
once the series is converged. The limiting factor here is 
the accuracy of the numerically obtained Chebychev ex- 
pansion coefficients. This effect becomes more severe, as 
the argument of the Chebychev polynomials, At/^E (and 
thus the largest expansion coefficient) becomes larger and 
larger. 

The errors obtained by the Chebychev propagator with 
iterative time ordering of the order of 10~^^ to 10~^ 
have to be compared to those obtained by the standard 
Chebychev propagator, i.e. neglecting all effects due to 
time ordering. The latter yields maximum solution er- 
rors, £^5''^, of the order of 10"'' for At = 10a.u.= T/900 
and 10-3 for At = 40 a.u. The smallest ef^'^ that can be 
achieved without time ordering is of the order of 10"^ for 
At = 10^2 a.u.= T/900000. Thus the numerical results 
obtained with the iterative method are highly accurate 
compared to those obtained by the standard Chebychev 
propagator neglecting time ordering. 

Regarding the numerical efficiency of the Chebychev 
propagator with iterative time ordering, several conclu- 
sions can be drawn from Table HI First of all, it is abso- 
lutely sufficient to choose the number of sampling points 
within the interval At, TVj, only slightly larger than the 
order of the expansion of the inhomogeneous term, m^.. 
Doubling Nt doesn't yield better accuracy but requires 
more CPU time. Second, we expect an optimum in terms 
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A -I- 


Nt 


mk 


A'^Chcby 




^ norm 


CPU time 


^max 


10 


6 


4 


10 


1.7 ■ 10"'^ 


1.1 ■ 10"" 


23 s 


3 




12 


4 


10 


1.7 ■ 10~'' 


1.1 • 10"" 


47 s 


3 


20 


7 


5 


11 


4.1 ■ 10"'' 


1.8 • 10"" 


14 s 


4 




14 


5 


11 


4.1 ■ 10~" 


1.8 • 10"" 


29 s 


4 


40 


7 


5 


14 


3.1 ■ 10"" 


1.2 • 10~" 


8s 


4 




14 


5 


14 


3.1 ■ 10~" 


1.2 • 10"" 


15 s 


4 


80 


8 


6 


16 


1.9 ■ 10"" 


1.1 • 10"" 


6s 


5 




16 


6 


16 


1.9 ■ 10"" 


1.1 • 10"" 


11 s 


5 


100 


9 


7 


17 


8.3 ■ 10"" 


4.0 • 10"" 


5s 


5 




18 


7 


17 


8.3 ■ 10"" 


4.0 ■ 10"" 


9s 


5 


300 


10 


8 


29 


1.9 ■ 10"'" 


1.0 • 10"'" 


3.4 s 


6 




20 


9 


29 


1.9 ■ 10"'" 


1.0 • 10"'" 


5.3s 


6 


600 


12 


10 


32 


5.7 ■ 10"'" 


3.1 • 10"'" 


2.6 s 


6 




24 


10 


32 


5.7 ■ 10~'" 


3.1 • 10"'" 


4.2 s 


6 


700 


12 


10 


33 


7.8 ■ 10"'" 


3.6 • 10"'" 


2.1s 


7 




24 


10 


33 


7.8 ■ 10 


3.6 • 10 


3.8s 


7 


800 


14 


12 


35 


5.2 ■ 10"'" 


2.3 • 10"'" 


2.5 s 


8 




28 


12 


35 


5.2 ■ 10"'" 


2.3 • 10 


4.3 s 


8 


900 


15 


13 


36 


1.1 • lO"'' 


5.3- 10"'" 


3.0 s 


8 




30 


13 


36 


1.1 • 10"^ 


5.3- 10"'" 


5.1s 


8 


1000 


17 


15 


38 


3.6-10"'' 


7.0 • 10"'" 


3.5 s 


9 




34 


15 


38 


3.6 • lO"'' 


7.0 • 10"'" 


5.8s 


9 


TABLE I. The maximum error of the solution, e^i^, and the 


maximum deviation of the norm 


from unity, e™omj occuring 



in the overall propagation time are hsted together with the 
required CPU time for several short time intervals At. Nt 
denotes the number of sampling points within At, A'^choby the 
largest number of Chebychev coefficients in the expansion of 
Fm, mfe the order of the expansion of the inhomogeneous term 
and fcniax the largest number of the iterations for time ordering 
occuring for all time intervals [t„,f„+i]. 

of CPU time as At is increased. A Chebychev expansion 
always comes with an offset and becomes more efficient 
as more terms in the expansion but less time steps are 
required (this concerns both Chebychev expansions, that 
for the inhomogeneous term of order uik and that for 
the time evolution operator, i.e. for the F™, of order 
A^chcby)- However, this trend is countered by a higher 
number of iterations for time ordering, fcmax- Accord- 
ing to Table [H the optimum in terms of CPU time is 
found for At « 700 a.u. Finally, the order required in the 
Chebychev expansion of the inhomogeneous term, m^, 
stays comparatively small, well below the values where 
the transformation between the Chebychev coefficients 
and the polynomial coefficients becomes numerically in- 
stable, cf. Appendix [K[ 

B. Driven harmonic oscillator 

As a second example, we consider a harmonic oscillator 
of mass m = 1 a.u. and frequency to = 1 a.u. driven by a 
linearly polarized field. The time-dependent Hamiltonian 
is given by 

H(r; t) = -^1^2 + + rEoSit) cosM , (35) 



where Eq is the maximum field amplitude, S{t) the shape 
function given by Eq. (|29p . and luq is the frequency of 
the driving field. The final time is set to T = 100 a.u. 
The Hamiltonian is represented on a Fourier gridi^ with 
A^grid = 128 grid points, and r^ax = 10a.u.= —r^m- The 
transition probabilities and expectation values of posi- 
tion and momentum as a function of time are known 
analytically fiii^ 

Taking the initial wave function \ip{t = 0)) to be the 
ground state of the harmonic oscillator, we again mea- 
sure the deviation of the numerical from the analytical 
solution, Esoi, and the deviation of the norm of {tpit)) 
from unity, Enorm, for the time-dependent probability of 
the oscillator to be in the ground state. The pulse am- 
plitude is chosen to completely deplete the population of 
the ground state. 

Two cases are analyzed which both correspond to 
strong resonant driving of the oscillator. In the first 
case the rotating-wave approximation is invoked, i.e. we 
set ojQ = 0. This eliminates the highly oscillatory term 
from the field, keeping only the time-dependence of the 
shape function (moderate time-dependence). In the sec- 
ond case, the rotating-wave approximation is avoided, 
luq = w, i.e. the time-dependence of the Hamiltonian 
is much stronger than in the first case (strong time- 
dependence). 

In order to compare the Chebychev propagators with 
and without time ordering, we first list the smallest e'^^i^ 
and elfomi achieved by the standard Chebychev prop- 
agator without time ordering in Table [Hi The stan- 
dard Chebychev propagator was developed for time- 
independent problems and is most efficient for large time 
steps. Here, however, extremely small time steps At have 
to be employed to minimize the error due to the time- 
dependence of the Hamiltonian. Consequently, the re- 
quired CPU times become quickly very large. Note that 
the deviation of the norm from unity is much smaller 
than the error of the solution. This means that norm 
conservation cannot serve as an indicator for the error 
due to the time-dependence of the Hamiltonian. This er- 
ror is clearly non-negligible even for the very small time 
steps shown in Table [TTl 

Tables IIIII and IIVI compare the results for the driven 
harmonic oscillator obtained by the Chebychev propa- 
gator with iterative time ordering (ITO) and without 
time ordering (standard Chebychev propagator). The 
rotating-wave approximation is invoked in Table IIIII 
cjQ = 0, and avoided in Table HVl ujq = u). 

In the case of the rotating-wave approximation, both 
propagators conserve the norm on the order of 10"^^. 
However, only the propagator with iterative time order- 
ing achieves an accuracy of the solution of the same order 
of magnitude while the standard propagator yields errors 
of the order of 10"^ for the time steps listed in Table Hill 
The smallest maximum error of the solution achieved by 
the standard Chebychev propagator for = is of the 
order of 10~® for a norm deviation of the order of 10"^^, 
cf. Tab. ini However, this requires a prohibitively small 
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At 




ni ax / \ 
enorm 1^0 = <-^j 


max / c\\ 
Enorm (^^0 = Oj 


essr ( = 


max/ r\\ 

eSf (^^0 = 0) 


CPU time 


10-" a.u. 
10-5 a.u. 
10-" a.u. 
10-^ a.u. 


4 
5 
7 
10 


6.8 • lO-'' 
6.6 ■ 10-1" 
1.3 ■ 10-" 
6.5 • 10-12 


6.7- 10-^ 

6.6 ■ 10-1" 

4.7 ■ 10-11 
7.3 ■ 10-12 


4.6 ■ 10-** 

4.7 ■ 10-^ 
4.7 ■ 10-" 
4.7 ■ 10-5 


6.7 ■ 10-" 
4.4 ■ 10-" 
4.2 ■ 10-* 
4.2 ■ 10-^ 


1 h 16 m 28 s 
9 m 26 s 
1 m 28 s 
12 s 



TABLE II. Driven harmonic oscillator with [uoo = 0) and without (ujo — to) the rotating-wave approximation for the standard 
Chebychev propagator without time ordering. A'^chcby is the number of Chebychev polynomials required for the expansion of 

g-iHAt 



with iterative time ordering (ITO) 


without time ordering (standard) 


At 


Nt 


rrik 


-^Chcby 


max 
^ norm 


max 


CPU time 


h 

"'max 


At 


A^Choby 


max 
^norm 


max 


CPU time 


0.01 a.u. 


10 


8 


9 


5.3 ■ 10-'-' 


5.3 ■ 10-1'' 


1 m 54 s 


2 


0.01 a.u. 


18 


1.4 • 10-12 


4.2 ■ 10-" 


2s 


0.02 a.u. 


10 


8 


10 


1.4 • 10-1^ 


1.4 • 10-1^ 


58 s 


2 


0.02 a.u. 


18 


5.7-10-12 


8.5 ■ 10-" 


1.28s 


0.04 a.u. 


10 


8 


15 


5.5 ■ 10-1^ 


4.5 • 10-1^ 


31s 


2 


0.04 a.u. 


29 


6.15 ■ 10-12 


1.7- 10-5 


Is 


0.1 a.u. 


10 


8 


24 


1.4 • 10-12 


1.4 • 10-12 


27 s 


3 


0.1 a.u. 


44 


1.5 ■ 10-12 


4.2 ■ 10-5 


0.52 s 



TABLE III. Comparison of the Chebychev propagator with iterative (ITO) and without time ordering for the driven harmonic 
oscillator in the rotating-wave approximation {uio — 0). Notation as in Table |T] 
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FIG. 2. (color online) Strongly time-dependent Hamiltonian 
(uo = <^): Comparison of the Chebychev propagators without 
time ordering (standard) and with iterative time ordering in 
terms of the difference between the numerical and analytical 
solution, esoi(t). 



time step, Ai = 10 a.u. 

Even for a very strongly time-dependent Hamiltonian, 
when the rotating-wave approximation is not invoked 
(luq ~ Lu), the Chebychev propagator with iterative time 
ordering yields similarly accurate results, with errors of 
the order of IQ-^^, cf. Table HVl For comparison, the 
error obtained for the standard Chebychev propagator 
without time ordering is of the order of 10"'^ for the time 
steps reported in Table \TV\ The smallest errors achieved 
with the standard Chebychev propagator are of the or- 
der of 10"^ for a norm deviation of the order of 10'^^ for 
extremely small time steps, cf. Table HIl 

The error of the solution, eso\{t), is shown in Fig. [2] 
as a function of time for different time steps and a very 
strongly time-dependent Hamiltonian, luq = lu. This il- 
lustrates the superiority of the Chebychev propagator 



with iterative time ordering in terms of accuracy. A 
comparison of Tables [TTl and lIVl reveals furthermore that 
the Chebychev propagator with iterative time ordering 
is also more efficient than a standard Chebychev propa- 
gator with very small time step if a high accuracy of the 
solution is desired. 

Since we have established the Chebychev propaga- 
tor with iterative time ordering as a highly accurate 
method for the solution of the TDSE with explicitly 
time-dependent Hamiltonian, it is worthwhile to com- 
pare it to alternative propagation methods for this class 
of problems. In the following we will consider the {t,t') 
method^ and a fourth-order Runge-Kutta scheme. The 
{t, t') method provides a numerically exact propagation 
scheme by translating the problem of time ordering into 
an additional degree of freedom of a time-independent 
Hamiltoniaufii The TDSE for the Hamiltonian in the ex- 
tended space is solved by numerically exact propagation 
schemes such as the Chebychev or Newton propagators<^ 
The {t, t') method is, however, relatively rarely used 
in the literature due to its numerical costs in terms of 
both CPU time and required storage. On the other 
hand, Runge-Kutta schemes are extremely popular in the 
literaturci^ They are potentially very accurate if a high- 
order variant is employed. Note that high order of a 
Runge-Kutta method implies evaluation of the Hamilto- 
nian at several points within the time step [t„, tn+i]- 

In order to achieve a fair comparison between the 
Chebychev propagator with iterative time ordering and 
the (t, t') method, first the parameters which yield an op- 
timal performance of the (t,t') method for our example 
have to be determined. The required CPU time and the 
errors, e™5"^ and efj^'^m: a function of the number of 
grid points, iVf and Nt, are listed in Table fVl In case of 
a moderate time-dependence of the Hamiltonian corre- 
sponding to the rotating-wave approximation (uj = 0), a 
fairly small number of grid points in both t and t' is suf- 
ficient. Note that the number of points in the auxiliary 
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with iterative time ordering (ITO) 


without time ordering (standard) 


At 


Nt 


nik 


-^Chcby 


max 
^ sol 


max 
^norm 


CPU time 


^max 


At 


-^Choby 


max 
^norm 
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^sol 


CPU time 


0.01 a.u. 


10 


8 


9 


8.2 ■ 10""" 


3.5 ■ 10-" 


2 m 5s 


3 


0.01 a.u. 


18 


1.5 • 10""" 


4.6 • 10-* 


2s 


0.02 a.u. 


10 


8 


10 


2.5 ■ 10-" 


3.6 ■ 10-"" 


1 m 34 s 


3 


0.02 a.u. 


23 


4.7 • 10""" 


9.3 • 10-* 


1.28 s 


0.04 a.u. 


10 


8 


15 


3.6 ■ 10"" 


5.5 ■ 10""" 


56 s 


3 


0.04 a.u. 


29 


8.9 • 10""" 


1.8 • 10"^ 
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TABLE IV. Comparison of the Chebychev propagator with iterative (ITO) and without time ordering for the driven harmonic 
oscillator without the rotating- wave approximation {ljo = lj). Notation as in Table U] 
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16 m 37s 
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13 
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13 


43 s 
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16 


906 


3.4 


10- 


13 


1.7 
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13 


32 s 
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8 


1740 


2.4 
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13 


1.2 


10" 


13 


30 s 
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43 


2.5 


10" 


lU 


3.2 


■ 10" 
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17 m 43 s 


2048 


2048 
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10" 


11 


2.8 
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11 


1 h 08 m 46 s 


2048 
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3.3 


10" 


11 


1.6 


10- 


11 


36 m 58 s 


2048 


256 


119 


3.3 


10" 


11 


1.6 


10- 


11 


30 m 36 s 


2048 


128 


204 


3.3 


10- 


11 


1.6 


10- 


11 


25 m42s 


2048 


64 


366 


3.3 


10" 


11 


1.6 


10" 


11 


22 m 58s 


2048 


32 


680 


3.3 


10" 


11 


1.6 


10" 


11 


21 m 21s 


2048 


16 


1295 


3.3 


10" 


11 


1.6 


10" 


11 


21 m 14s 



TABLE V. Performance of the {t, t') method. The required 
CPU time and the errors, elfol-m a-nd e^j"^, are listed for dif- 
ferent numbers of sampling points of the t' coordinate, Nti 
and different numbers of sampling points within [0, T], Nt to- 
gether with the number of required terms in the Chebychev 
expansion, Achoby. The upper (lower) part corresponds to 

LJO = {UJO = Uj). 



■■ 1 1 1 1 1 1 1 1 1 ' 


"1 1 ' ' ! 


r ♦— ♦ Chebychev with iterative time ordering \ * 




: •-• (t,t')-method N 




i-^ Chebychev without time ordering 




; »— • Runge-Kutta of 4th order 

" 1 1 1 1 1 1 1 1 1 . ... 


..1 . . .." 



Ixl0"'^xl0"''lxl0 "lxlO"'°lxlo '' IxIO'" 1x10"' 1x10"' 1x10"' 1x10"^ 1x10" 



FIG. 3. (color online) Comparison of propagation methods 
for strongly time-dependent Hamiltonian (ljo = <^) in terms 
of the CPU time required in order not to exceed a given max- 
imum error of the solution, e^^^. 
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^sol 


max 
^norm 
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ITO 
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RK4 
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1.2 • 10-"^ 
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31s 
30 s 
38 m 24 s 




ITO 
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1.7- 10-"^ 
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21 m 14s 
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TABLE VI. Comparison of highly accurate methods. 



coordinate, Nti , is not known a priori. 

For a strong timc-dcpcndcnce, i.e. = a fairly 
large number of points for the auxiliary coordinate, t', is 
required, Nti = 2048. However, since the actual prop- 
agation involves a time-independent Hamiltonian, large 
time steps can be taken for the Chebychev propagator, 
resulting in the most efficient solution when Nt is small, 
Nt ~ 16, and correspondingly the number of Chebychev 
terms, Ncheby, is large. 

Table IVII reports the comparison between the Cheby- 
chev propagator with iterative time ordering (ITO), the 
{t, t') method, and the fourth-order Rungc-Kutta scheme 
(RK4) for the resonantly driven harmonic oscillator. For 
a moderate time-dependence, i.e. in the case of the 
rotating-wave approximation {ujq = 0), the {t^t') method 
and the Chebychev propagator with iterative time order- 
ing yield a similarly good performance in terms of er- 
rors and CPU time. Contradicting the common percep- 



tion of the Runge-Kutta scheme as a particularly efficient 
method; the CPU time for our example is found to be al- 
most two orders of magnitude and the error three orders 
of magnitude larger than for the Chebychev propagator 
with iterative time ordering and the (t, t') method. 

For strong time-dependence, i.e. resonant driving 
without the rotating-wave approximation [loq = uj), the 
Chebychev propagator with iterative time ordering is 
found by far superior in terms of both efficiency and accu- 
racy compared to the (t, t') method and the fourth-order 
Runge-Kutta scheme. In both cases, the Runge-Kutta 
scheme is the least accurate method. The smallest error 
of the solution, e™^'^, achieved with RK4 is of the order 
of 10"^ for ajQ = oj and At = 10"^ a.u. and of the order 
of IQ-^ for cjQ = with the same At. We have not tested 
smaller time steps, since already with At = IQ-^a.u., 
RK4 is the least efficient of the three methods in terms 
of CPU time. 

Figure [3] illustrates how much CPU time is required 
for a given maximum error of the solution, s^^^. A clear 
separation between highly accurate methods (Chebychev 
propagator with iterative time ordering, (t, t') method) 
and less accurate methods (standard Chebychev propa- 
gator without time ordering, fourth-order Runge-Kutta 
scheme) emerges. If a highly accurate method is desired, 
the Chebychev propagator with iterative time ordering 
appears to be the best choice. It outperforms the (t,t') 
method not only in terms of CPU time as shown in Fig. [3] 
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but also in terms of required memory. In the intermedi- 
ate regime realizing a comprise between accuracy and 
efficiency, the Chebychev propagator with iterative time 
ordering is still the best choice. While the less accurate 
methods that ignore time ordering become prohibitively 
expensive, the {t, t') method does not cover this regime. 
This is due to the choice of Nt' - if it is large enough, the 
calculation is converged and the error is very small, if it 
is too large, convergence cannot be achieved and norm 
conservation is violated. Only for cases, where a lim- 
ited accuracy of the solution is sufficient {Ssoi > 10^^), 
the standard Chebychev propagator and the fourth-order 
Runge-Kutta scheme represent the most efficient propa- 
gation schemes. 



C. Wave packet interferometry 

Our third example applies the Chebychev propagator 
with iterative time ordering to a model that cannot be in- 
tegrated analytically. It explores the effect of time order- 
ing on phase sensitivity as employed in coherent control. 
Wave packet interferometry has first been demonstrated 
in the early 1990si^ A pair of electronic or vibrational 
wave packets are made to interfere by two laser pulses. 
This represents a conceptually very simple prototype of 
quantum controli^ The interference is controlled by the 
relative phase between the two pulses. 

Our example is inspired by a recent experiment We 
consider two harmonic oscillators that are coupled by a 
laser field. 



H 



t + Ve(r) 



(36) 



where T denotes the kinetic energy and 

1 



Ve(r) 



2m 
1 

2m 



For simplicity we again take m = 1, = We = 1, and 
/i = 1 a.u. The Hamiltonian is represented on a Fourier 
grid with A^grid = 128, Tmin = —10 a.u., r,„ax = 12 a.u. 
and Te = 3.5 a.u. Starting from the vibronic ground 
state, a pump pulse is applied to create a wave packet 
in the excited state, cf. Fig. U^. The excited state wave 
packet oscillates back and forth in the excited state po- 
tential with a period of 2tt /uje- The control pulse, with 
parameters identical to those of the pump pulse, can be 
applied with different time delays. If it is applied after 
one vibrational period, a relative phase equal to zero in- 
duces constructive interference while a relative phase of tt 
induces destructive interferencei^ Different time delays 
combined with a different choice of the relative phase 
yield the same result.— Constructive interference implies 
an increase of population in the excited state, while for 
destructive interference the wave packet is deexcited to 
the ground state. 




Pi 



1 ' 1 ' 1 ' 1 ' 


1 1 1 1 1 1 _i 


y constructive \ 




Y interference \ 


S. destructive 




\ interference 


T . 1 . 1 . 1 , 


1 .^"T"^. 1 



-90 -45 45 90 135 180 225 

relative phase (p [degrees] 

FIG. 4. (a) Schematic representation of the generation of 
wave packets in the excited state with an ultrashort laser 
pulse, (b) Ratio of the population on the excited state as 
a function of the relative phase between the pump and the 
control pulse, ti is the end of the pump pulse. 



The excited state population that was measured in the 
experiment by a probe pulse^ can be simply calculated, 
KV'elV'e)!'^- The ratio of excited state population at the 
final time T and at time ti, just after the pump pulse. 



i?(^) 



|(^e(r)|^e(r))p 
\{MtiMe{h)W 



(37) 



depends on the relative phase between the two pulses, Lp. 
This dependence is illustrated in Fig. ^p. For Lp = Q. the 
population increases by a factor of fotir and for = tt 
complete de-excitation is observed. 

Since an analytical solution is not available for this 
example, we take the solution obtained by the Chebychev 
propagator with iterative time ordering as the reference. 
The accuracy of propagators without time ordering is 
analyzed in terms of the relative error e^jjj, 



\RiToiv) - Riip)\ 
RiToiv) 



(38) 



They are shown for the standard Chebychev propagator, 
the split propagator and the fourth-order Runge-Kutta 
scheme in Figs. [5] and [6] for different pulse energies (re- 
spectively, pulse areas) and At = 10~^ a.u. (which has to 
be compared to the duration of the pulse, 0.3 a.u. and the 
vibrational period, 27ra.u.). Fig. [5] corresponds to (al- 
most) constructive interference, (/s « 0, Fig.[6]to (almost) 
destructive interference, « tt. Overall, the relative 
errors obtained are smaller for wave packet interference 
compared to the examples of the previous sections IIV Al 
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FIG. 5. (color online) Constructive wave packet interfer- 
ence: Accuracy of the standard Chebychev propagator, the 
split propagator and the 4th order Runge-Kutta scheme with 
respect to the Chebychev propagator with iterative time or- 
dering for different pulse areas. 



timc-ordcring and the split propagator appear to be only 
weakly affected by time ordering for constructive intef- 
erence. The results obtained with these two propagators 
for destructive interference are much more sensitive to 
time ordering effects and relative errors reach between 
10"^ and lO^"' for weak and strong pulses, respectively. 
The error of the split propagator is due to two effects 
- time-ordering and the non-vanishing commutator be- 
tween kinetic and potential energy while the error of the 
standard Chebychev propagator is solely due to time or- 
dering. For weak pulses, the standard Chebychev prop- 
agator yields more accurate results than the split prop- 
agator, cf. Fig. [6l However, for strong pulses (pulse 
area of tt/2) roughly the same accuracy is achieved by 
the standard Chebychev propagator and the split prop- 
agator. This indicates that the neglected time-ordering 
becomes the dominating source of error. 



CONCLUSIONS 
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FIG. 6. (color online) Destructive wave packet interference: 
Accuracy of the standard Chebychev propagator, the split 
propagator and the 4th order Runge-Kutta scheme with re- 
spect to the Chebychev propagator with iterative time order- 
ing for different pulse areas. 



and IIVBI We attribute this to the fact that the pump 
and control pulse arc very short compared to the vibra- 
tional time scale of the oscillators. In this regime of im- 
pulsive excitation, the pulses act almost as (S-functions, 
and there is not enough time to accumulate large errors 
due to neglected time ordering. However, even in this 
regime the errors are non-negligible. As expected the er- 
rors become larger with increasing pulse intensity. The 
Runge-Kutta scheme yields similar errors for both con- 
structive and destructive interference. The results ob- 
tained by the standard Chebychev propagator without 



We have developed a Chebychev propagator based on 
iterative time ordering to solve TDSEs with an explic- 
itly time-dependent Hamiltonian. The key idea consists 
in rewriting the term of the TDSE that contains the 
time-dependence of the Hamiltonian as an inhomogene- 
ity. This inhomogeneity can be approximated iteratively. 
At each step of the iteration, the Chebychev propaga- 
tor for inhomogencous Schrodingcr equations^ is em- 
ployed. Convergence is reached when the wave functions 
of two consecutive iteration steps differ by less than a 
pre-specified error. Time ordering is thus accomplished 
in an implicit manner. 

We have outlined the implementation of the algorithm 
and demonstrated the accuracy and efhciency of this 
propagator for three different examples. A comparison 
to analytical solutions and other available propagators 
has shown our approach to be extremely accurate, yet 
efficient, in particular for very strong time dependencies. 

The importance of correctly accounting for time or- 
dering effects^i has been demonstrated for destructive 
quantum interference phenomena. In most of the litera- 
ture on coherent control, accurate propagation methods 
are employed but time ordering effects are completely 
neglected. This is not justified, in particular for appli- 
cations such as optimal control theory or high-harmonic 
generation where the fields are very strong. 

The approach of rewriting parts of the TDSE as an in- 
homogeneity can be extended to other classes of problems 
where numerical integration is difficult. An obvious ex- 
ample is given by non-linear Schrodinger equations such 
as the Gross-Pitaevski equation. By rewriting the non- 
linear term as an inhomogeneity. it should be possible to 
derive a very stable propagation scheme. 
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Appendix A: Transformation to obtain the coefFicients 
|$(^ )) from the Chebychev expansion coefFicients $j) of 
the inhomogeneous term 



In order to make use of Eq. (|14p . a transformation link- 
ing the Chebychev coefficients, \<tj) that are calculated 
by cosine transformation of the inhomogeneous term, cf. 
Eq. (fT3|) . to the coefficients |<I>*^-' ^) appearing in the for- 
mal solution of the inhomogeneous Schrodinger equation, 
cf. Eqs. (fT5|) - PT|) . is required. Assuming r e [0,t], then 
f = 2T/t — 1, and Eq. p4)) becomes 



rn—l 



3=0 



j'=0 



3' 



(Al) 



Replacing t by r in the right-hand side of (|A1[) . one ob- 
tains 



m — 1 



3=0 



3'=0 



(A2) 



The Chebychev polynomials can be expanded in powers 
of f , 



P3ir)=EC3/j^ 



k=0 



(A3) 



P,-+i(f) = 2fP,(f)-P,_i(f), 



(A4) 



the coefficients Cj_k satisfy a corresponding recursion re- 
lation, cf. Eq. (A6) of Ref. [H. Inserting Eq. into 
Eq. ([XT]) yields 



m-l 3 



m-l j 



j=0 k=0 j'=0 k=0 

Introducing 



(A5) 



c. 



3-k, 



k\{j' - k)\ 21 



\0^3,k) = 
\p3'..k) = 

Eq. (|A5|) is rewritten 

m — 1 J m — 1 J 

E E 1".^)^'= = E E \Pr^-)^' 



(A6) 



Calculation of the ') from the is thus equivalent 
to calculate the \Pji.k) from the \oij^k)- Note that the 
powers of f in Eq. (|A6P occur in the inner sums. In order 
to transform them to the outer sums, first the left-hand 



m-l 3 

EE 

j=0 fc=0 
m-l J 



„ 1 :i ?ri— i 

^i$,)f'= = |ao,o> + E l"!.'^-)^' + E l"2,fe>^' + • ■ ■ + E l«™-i.fc>^* 



A:! 



EE%^|^.) 



j=0 fc=0 

m-l 3 

EE^i'J'^-^^ 

j=0 fc=0 



m-l ?n— 1 



fc=0 



fc=0 



m — 1 



E + E + E l"i-2)^^ + • • • + |a,„_i,™_i)T™ 1 , 



3=0 j = i 

1 m— 1 

E E i".-.^)^^' 



J=2 



(A7) 
(A8) 
(A9) 



Similarly, the right-hand side of Eq. (|A5[1 is written 



m— 1 m— 1 



777. — Xj^ ^ ''^ — -L ' ' — -L 



j'=0 k=0 ' j'=0 k=j' 



(AlO) 



Equating the right-hand sides of Eq. (fA9|) and Eq. (jAlOp . 

the \(}j',k) are obtained, 

m — 1 m-l 

Y.\f^r,k)^Y.\^3,k), 0<k<m-l. (All) 

j' = fe j=k 
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Replacing \aj^k) and by their definition yields 



m— 1 



^ k\{j' -k)\ 23 
j —k 



^ 7n—l 



1 P 



m-l j> 



(j'-fc)!2^' 



This leads to the hierarchy of equations 



rn — l 



j-m—l 



m-l „ 

E^l^i). 0<fc<m-l, 

j=k 

7n — 1 

EQ,fc|$j>, 0<fc<m-l. 
j=fe 



)m — 2 



(A12) 



(A13) 

j'=k j=k 

The coefficients |$'^ •*) can thus be determined step by step from the coefficients of the Chebychev expansion, 



om — 1 
|^(m-l)N ^ 



m— l.m — 1 l^m— 1 



r 



fc = m - 2,0. 



(A14) 
(A15) 



Note that the transformation given by Eqs. (jAl4p and 
(jA15|) becomes numerically instable for large orders, 
TO ~ 100. In our applications, such a large to would cor- 
respond to time steps larger than the overall propagation 
time and was never required. 
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